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ABSTRACT 

Using magnet ohydro dynamic (MHD) adaptive mesh refinement simulations, we study the formation 
and early evolution of disk galaxies with a magnetized interstellar medium. For a 10^^ M© halo with 
initial NFW dark matter and gas profiles, we impose a uniform 10~^ G magnetic field and follow its 
collapse, disk formation and evolution up to 1 Gyr. Comparing to a purely hydrodynamic simulation 
with the same initial condition, we find that a protogalactic field of this strength does not significantly 
influence the global disk properties. At the same time, the initial magnetic fields are quickly amplified 
by the differentially rotating turbulent disk. After the initial rapid amplification lasting ~ 500 Myr, 
subsequent field amplification appears self-regulated. As a result, highly magnetized material begin 
to form above and below the disk. Interestingly, the field strengths in the self-regulated regime agrees 
well with the observed fields in the Milky Way galaxy both in the warm and the cold HI phase and 
do not change appreciably with time. Most of the cold phase shows a dispersion of order ten in 
the magnetic field strength. The global azimuthal magnetic fields reverse at different radii and the 
amplitude declines as a function of radius of the disk. By comparing the estimated star formation 
rate (SFR) in hydrodynamic and MHD simulations, we find that after the magnetic field strength 
saturates, magnetic forces provide further support in the cold gas and lead to a decline of the SFR. 
Subject headings: galaxies: formation — galaxies: ISM 
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1. INTRODUCTION 

Rapid developments in observational cosmology re- 
cently established a stan dard mode l of cosmology, the 
ACDM model (Spergel et al. I l2007h . In principle, this 
makes the problem of galaxy formation essentially an 
initial value problem. Indeed, Gpc-scale N-body simu- 
lations of ACDM cosmology are revealing the details of 
large-scale structure formation and clustering;, matching 
observations in many aspects (ISpringel et al. I [2005). In 
this framework, it should be possible to understand the 
formation and evolution processes of individual galaxies. 
For this purpose, it is crucial to understand the structure 
of interstellar medium (ISM) and the dominant physics 
controlling star formation. 

Analytical and semi-analytical studies of disk galaxy 
for mation are useful in understanding the global proper- 
ties (iMo et al.llT99 8: Efstathiou 2000; Dutt on et al.ll2007l : 
iKampako glou Silk i2007i : tStringer & Benson] |200^. 
Because of the complex nature of galaxy formation and 
the wide range of physical processes involved, numerical 
simulations also prove particularly useful. Hydrodynam- 
ical simulations of galaxy formation in both cosmological 
and isolated setups have taught us a lot about the de- 
tailed structure of the galactic ISM and star formation 
(see e.g. Kravtsov 2003, Li et al. 2005, Springel & Hern- 
quist 2005, Okamoto et al. 2005, Tasker & Bryan 2006, 
Tassis et al. 2006, Governato et al. 2007, Wada & Nor- 
man 2007, Kaufmann et al. 2007, Robertson & Kravtsov 
2007 for recent studies). However, there are still many 
physical processes we need to include to make the simu- 
lations realistic. 
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In this work, we focus on magnetic fields which have 
not yet been included in previous galaxy formation sim- 
ulations. Magnetic fields may significantly influence the 
structure and evolution of ISM and has been studied 
extensively previously by local magnet ohydro dynamic 
(MHD) simul ations of the ISM (Mac Low et al. 20^; 
Balsara et al.l l2004t Ide Av illez & Breitschwerdt 2005; 
Piontek fc Qstriker II2007I : iHennebelle Inutsuka ..2006> ) 
and by observations (e.g. Beck 2007, Crutcher 1999). 
On a larger scale, the effect of magnetic fields is less 
clear because co smological MHD sirn ulations are still at 
an early stage. iMiniati et al.l (|200ll ) performed cosmo- 
logical hydrodynamic simulations for a passive magnetic 
field. Cosmological MHD simulations have been done in 
SPH (e.g. Dolag 1999) in the context of galaxy cluster 
formation and Eulerian-code simulations are also begin- 
ning to be explored (e.g. Li et al. 2006). None of those 
calculations have as yet resolved the internal structures 
of disk galaxies. 

Besides the potential importance in the dynam- 
ics of galaxy formation, galactic magnetic fields are 
a key ingredient in many other astrophysical prob- 
lems. For example, they play a domina nt role i n the 
origin and propag ation of cosmic rays (|Fermi I Il949l : 
IStrong et al.l l200"7n . which may also be dynamic ally im- 
port ant in galaxy formation (Fermi 1954; EnB lin et al.l 
l2007l ). Furthermore, magnetic fields may b e important 
in th e dynamics of supernova remnants (Balsara et al.l 
[200Th . in regulating the internal tur bulence of giant 
molecular clou ds by dr iving outflows (|Li fc Nakamural 
i2006; Banerje e fc Pudri tz 2007) and in the colla pse of 
individual molecular cloud co res ( Mouschovia sI 1 19871: 
Balsara et aDl200ll : iTillev fc Pudritz 2007: Mellon fc Li I 



20071 ). 
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While the importance on the ISM is generahy ac- 
cepted, the origin of the galactic magnetic field is still 
an unsolved problem (see e.g. Kulsrud & Zweibel 
2007 for a recent review). The traditional theory 
of a galactic dynamo, the alpha- Omega dynamo, is 
based on the mean fiel d dynamo theory introd uced by 



ISteenbeck et al. I (119661) and later dev e loped bv iParker 
(|1971[) and IVainshtein fc RuzmaikiiTI (|1971f ). H 



iFerriere 



(|1992l ) has extended this galactic dynamo theory to 
include the expansion of supernova bubbles. Tur- 
bulent a mplification in protogalaxies have also bee n 
studied (|Pudritz Silk I IT^M iKulsrud et aT~l 119971 ). 
Those theories are built on kinematic dynamo. A 
different paradigm for disk dynamo comes with the 
discovery of the magneto-rotational instability (MRI) 
(Chandrasekhar l [T9p; [Balbus & Hawley ^ HImI). In 
MRI, field amplification is not only a consequence of the 
turbulent velocity field, it also produces the turbulent ve- 
locity field itself. It is still unclear what is the dominant 
mechanism that amplifies the galactic magnetic field. In 
contrast to the case of dynamos in planets , because of the 
enormous inductance of the ISM (iFermi 1 11949). there is 
no decaying problem for the galactic field. So the essen- 
tial part of the problem here is indeed only the initial 
amplification. 

In the following, we first describe the numerical al- 
gorithms and initial conditions. Then in section 3 we 
present the results. Sections 4 and 5 are devoted to dis- 
cussion and conclusions, respectively. 

2. NUMERICAL METHODS AND MODELS 

2.1. MHD Algorithm 

We have developed a 3D adaptive mesh refinement 
(AMR) MHD code based o n the cosmological AMR 
hydrodynamics code Enzo (|Bryan fc NormaiTI [l997i : 
Ib'Shea et al. 2004). We use the Dedner formulation 
of MHD equations (jDedner et al. 1120021 ). This conser- 
vative formulation of the MHD equations uses hyper- 
bolic divergence cleaning. It has als o been used in sev- 
eral other recent AMR MHD codes (Matsumoto "2006; 
Anderson et al. 2006). The conservative system is dis- 
cretized by method of lines. The Riemann solver and 
reconstruction sche mes have been tested extensively in 
I Wang et al. I (|2007l ) for relativistic systems. Since our 
Riemann solvers and reconstruction schemes are de- 
signed to work for any conservative systems, when ap- 
plying it to the Dedner MHD equations, we only redefine 
the conservative variables. Using these high-resolution- 
shock-capturing (HRSC) schemes designed for conserva- 
tive systems ensures the conservation of mass, momen- 
tum and energy, as well as obtaining correct shock posi- 
tions and velocities (|LeVeque Il2002f ). Interested readers 
can find mo re detailed descripti ons of our numerical al- 
gorithms in I Wang et al. I (|2007l ) and references therein. 
In this work, reconstru ction is do ne using the piecewise 
linear method (PLM) (jVan Leer I ' 1979). Fluxes at ceh 
interfaces are calculated using the local Lax-Friedrichs 
Riemann solver (LLF, Kurganov & Tadmor 2000), which 
is free of the carbuncle artifacts pr esent in soni e other 
contact-capturing Riemann solvers (Q uirk I [19941 ). Time 
integration is performed using the total variation dimin- 
ishing; (TVD) second order Runge-Kutta (RK) scheme 
(|Shu fc Osher I [198 8). The code inherits all of Enzo's 
parallel AMR capability as well as self-gravity, chem- 



istry, cooling and radiative transfer routines. We have 
tested the code using various standard hydrodynamics 
and MHD test problems including the ID Brio-Wu prob- 
lem, 2D MHD rotor, 2D Hydro and MHD Rayleigh- 
Taylor problem, 3D Sedov- Taylor problem, 3D Larson- 
Penston isothermal collapse problem, etc. 

2.2. Cooling 

A cooling function is used to calculate the radiative 
energy losses down to temperature 300 K. We use the 
cooling function fitted by Sarazin White I (|l987r ) for 
T > 10"^ K and bv lR.osen fc Bregman I (|l995f ) for 300 K < 
T < 10"^ K. 

In global simulations of galactic disks, current compu- 
tational power does not allow to resolve the Jeans length 
down to molecular core scales where star formation ac- 
tually occurs. This may gives rise to artifi cial frag:menta- 
tion if one underresolves the Jeans length ([Truelove et al.l 
Il997f ). This potentially worry is avoided using a density 
dependent temperature fioor ensuring that a cell always 
resolves half of the local Jeans length (e.g. Wada & Nor- 
man 2007): 

284 . J ^Xk. (1) 
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Another common approach would be to create a colli- 
sionless star particle if a cell fails to resolve Jeans length 
(e.g. Tasker & Bryan 2006). 

2.3. Initial Conditions 

From pure N-body simulations it has been found that 
cold dark matter halos hav e a universal density profile 
well fitted by the NFW form (|Navarro et al. II1996D . Here 
we consider a small halo formed at redshift 2, with a 
mass of Myir = 10^^ M0and correspondingly, a virial 
radius Ryir = 21.5 kpc. The spin parameter is assumed 
to be A = 0.05 and concentration is Cyir = 10. Such a 
small halo allows for higher spatial resolution. By us- 
ing an isolated model, we are implicitly assuming that 
the effect of major and minor merger on the growth 
of galaxy disk is subdominant, which may be a reason- 
able assum ption if the halo is smaller than the Milky 
Way halo (|Guo White I l2007f ). However, our model 
galaxy is still larger than the small galaxies with rota- 
tion velocities < 30 km/s in which star formation may be 
strongly suppressed by the intergalactic UV background 
(jThoul fc WeinberdI 19961 ). The gravitational interaction 
between dark matter and baryons is a secondary effect. 
Consequently in this work we model the dark matter halo 
as a static external potential. 

The initial gas density profile follows also the NFW 
profile with its amplitude determined from the assumed 
baryon mass fraction = 0.1. The initial gas temper- 
ature is set by solving for local hydrostatic equilibrium. 
The rotation velocity is se t to be the same as that of 
the dark matter following S pringel fc Whitel (|T999). In 
addition, we also add to the gas random velocities of the 
same order as the dark matter virial velocity to model 
crudely the turbulent motions in hierarchical formation 
of protogalaxies. We put the NFW halo at the center 
of the simulation box with virial radius 0.1 (21 kpc) of 
the simulation box length. This makes sure that the gas 
evolution inside halos will be minimally affected by nu- 
merical boundary effects. The gas outside the halo virial 
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Fig. 1. — Face-on and edge-on density- weighted projections of density (left), temperature (middle) and thermal pressure (right) at 
t = 1.088 Gyr for the Hydro (top two rows) and MHD (bottom two rows) runs. The physical length of the plot is 11 kpc. 



radius is set to be the cosmic mean value at redshift 2 
with a temperature of T = 10^ K. The topgrid resolution 
is 64^, and there are three static nested grids refined by 
a factor of two around the halo. Then we let the code 
dynamically refine to higher levels according to the Jeans 
criterion with Jeans number 4. After the disk forms, al- 
most all the disk is refined to this highest level because 
it is strongly self-gravitating. The simulations were run 
with a maximum refinement level of six and seven, cor- 
responding to resolutions of 52 and 26 pc, respectively. 
We did not find noticeable differences in the results. So 
in the following, only the results of the 26 pc resolution 
run are discussed. 

The initial magnitude and topology of the magnetic 
field is uncertain. Cosmological simulations for the gen- 
eration of magnetic seed field will be required to address 
this important question. A tiny seed magnetic field of 
the order IQ"^-*^ G is always^ua ranteed by the Biermann 
battery effect ([Biermann I [l95Qi) . However, there are 



both observational and theoretical arguments suggest- 
ing larger pregalactic field strengths. For example, Fara- 
day rotation has been detected in high reds hift damped 
Lyman alpha system (jOren Wolfe iri995h . Also, the 
abundance of beryllium and boron in old Galactic halo 
stars is directly proportional to their iron abundance. 
Big bang nucleosynthesis does not produce beryllium nor 
boron. Hence, they may dominately be produced by spal- 
lation of low energy (tens of MeV) carbon and oxygen 
cosmic rays, suggesting that magnetic fi elds and cosmi c 
rays are already present at early times (Zweibel I I2QQ3I ). 
On a more theoretical ground, in hierarchical structure 
formation, any halo formed at late times must have had 
progenitors that hosted prior generations of stars. So 
any magnetic field produced by those stars, their super- 
novae, or their pulsar remnants would be amplified by the 
turbulent ISM in the galaxies containing them. Turbu- 
l ence driven by mergers may also amplify ma gnetic fields 
(jPudritz fc Silk |[T989l : [K^lsrud fc Anderson 1119921 Fur- 
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Fig. 2. — Velocity and Magnetic field plots at t = 1.088 Gyr. From left to right: 1. velocity field vectors overplotted on a density slice; 
2. magnetic field vectors overplotted on a density slice; 3. slice of plasma beta; 4. slice of . All face-on slices are through z = 0.5 and 
edge-on slices are through y = 0.5. The physical length of the plot corresponds to 11 kpc. 



thermore, iRees I (|2QQ6l ) argued that the magnetic fields 
from supernova eject a and extended radio lobes may 
build a field in excess of 10~^ G seed field for galactic 
disk forming in the late Universe. In this work, we adopt 
the simple assumptions that the initial magnetic field is 
weak, uniform and directed along the rotation axis. We 
will discuss the results of two simulations, with B = 
(Hydro run) and B = 10~^ G (MHD run) throughout. 
Initially, the magnetic field is dynamically unimportant, 
it will be twisted and amplified by the turbulent velocity 
field. 

This first study of MHD models of disk galaxy forma- 
tion only include a minimal set of physical processes rel- 
evant for the rich interplay between magnetic fields, the 
dynamics of disk formation and a multiphase ISM. We 
do not yet include supernova feedback, radiative heat- 
ing, cosmic rays, molecular cloud chemistry, ambipolar 
diffusion nor other physical processes that may be rele- 
vant. With those limitations in mind, our model should 
be viewed as a numerical experiment rather than a re- 
alistic simulation of current day disk galaxies. Some of 
those neglected feedback processes, e.g., supernova feed- 
back, may significantly change the conclusions drawing 
from the current calculations. 



3. RESULTS 

Fig. [1] shows the face-on and edge-on projections of 
the density, temperature and the thermal pressure. Fig. 
[2] shows face-on and edge on slices of velocity field and 
magnetic field lines overplotted on density, plasma beta 
f3 = p/ {B^ / (Stt)) and B^. Both figures correspond tot = 
1.088 Gyr. In some of the plots below, we will consider 
only the gas in the disk defined by a density threshold 
criterion p > 1 x 10~^^ g cm~^. 



3.1. Disk Structures 

As can be seen in Fig. (TJ the global disk structures are 
qualitatively similar in both the Hydro and MHD cases. 
The density projections show filamentary structures and 
high density "blobs" located along those filaments. Were 
we to Include molecular hydrogen chemistry and cooling, 
these blobs would presumably resemble giant molecular 
cloud complexes. This is similar to the density distri- 
butions of HI gas an d GMCs in obs ervations of nearby 
galaxies such as M33 (|Engargiola et al. 2003), which also 
show good correspondence between the filaments and the 
locations of the GMCs. The temperature projections 
show that those filaments are all cold, with T ^ 300 K. 
The lower density disk material has a temperature of 
T ^ 10^ K. In the MHD run. Fig. [2] shows that much of 
the cold gas is already strongly magnetized. 

In a self-gravitating disk, the vortex mode is 
the fastest unstable mode with growth rates 2 — 
3 times larger than the spiral den sity wave mode 
(|Mamatsashvili fc Chagelishvili I [20071 ). In an isolated 
non-rotating fiuid, the perturbation mode is divided into 
vortical and divergent types according to the Helmholz 
decomposition. And it is mainly the energy in divergent 
fiows that will be dissipated in shocks. In a differentially 
rotating disk, the vortex mode is a combination of vor- 
tical and divergent motions. Thus the morphology and 
statistical properties of large scale shocks which form the 
filaments are different in those two cases. However, there 
seems to be no quantitative comparative studies on this 
problem that we are aware of. 

To see how self-gravity affects the growth of vor- 
tex modes, we have carried out experiments with only 
the external NFW potential, neglecting the self-gravity 
of the gas. Analytical calculations in the thin sheet 
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Fig. 3. — Disk profiles for the Hydro run (thin fines) and the MHD run (thick fines) at t = 1.088. Top left gives the mass profile, top right 
the surface density profile, bottom left the velocity dispersion, bottom right the rotational velocity profile. In the mass profile plot, the 
thin dashed line is the mass profile of the static NFW h alo. In the surface density p rofile plot, the dashed line is the threshold gas surface 
density for star formation 5 M0pc~^ as observed bv lMartin &: Kennicutt"! (|200in . In the rotational velocity plot, the thin dashed line 
is the rotation velocity corresponding to the static NFW halo Vfj^^Npwir) = [^^ArFH^(^)/''^]"'^^^- 
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Fig. 4. — Toomre Q profiles for Hydro run (thin line) and MHD 
run (thick fine) cut t = 1.088. 

approximation predict that the growth rate of vortex 
mode is orders of magnitude smaller w ithout self-gravity 
(|Mamatsashvili fc Chagelishvili I l2QQ7l ) . Indeed in this 
case we saw that the disk shows much less filamentary 
structures and instead develops tightly wound spiral den- 
sity wave patterns. Thus the vortical turbulent motion 
in the fiducial cases is clearly driven by self-gravity. An- 
other interesting implication of this experiment is that 



many local galaxies do show tightly wound spiral pat- 
terns. Hence our experiment suggests that the gravita- 
tional force in the gas disks of those galaxies is domi- 
nated by the stars. This implies that for disk galaxies 
at high redshift that are still self-gravitating, the vor- 
tex modes may be the dominant structures regulating 
molecular cloud and star formation within them. In the 
following discussion we will call the turbulent velocity 
field seen in our simulations gravity-driven turbulence 
(|Gammie 1 l2(M IWadi^ et al. ll2QQ2h . 

As also can be seen in Fig.[Tl the MHD run has a some- 
what thicker disk than the Hydro run. This shows at 
t = 1.088 Gry, magnetic pressure is already large enough 
to provide additional support to the vertical force bal- 
ance, which can also be seen from Fig. [2l 

The pressure projections, show that the cold fila- 
ments have higher pressure than the warm gas and the 
cores of the filaments have even higher pressure. The 
overpressure in the cores is expected as they are self- 
gravitating (|LarsQnlll98lf ). The higher pressure in non- 
selfgravitating cold gas is consis tent with previous simu - 
lations of isolated disk galaxies (|Tasker fc Brvan1l2QQ6l ). 
This result is different from the idea that cold gas is 
formed by is obaric thermal instability in the galactic disk 
([Field II1965 V One important difference with the original 
thermal instability analysis is that instead of cosmic ray 



6 



heating, the heating sources in our simulation are pdV 
work, numerical viscous heating and numerical resistivity 
in the MHD run, all controlled by the gravity-driven tur- 
bulence velocity field in the disk. Thus in our simulations 
show that gravity-driven turbulent heating and radiative 
cooling will create a non-isobaric two phase medium. In 
the old isobaric theory, there is no prediction for the 
morphology and fraction of the cold phase. But in the 
turbulent heating scenario, one would expect the cold 
gas to have filamentary structures, which are naturally 
produced by a supersonic turbulent velocity field. 

Fig. [3] shows the azimuthally averaged radial profiles 
of mass, surface density, velocity dispersion and rota- 
tional velocity for both simulations at t = 1.088 Gyr. 
The Hydro and MHD calculations have similar density 
profiles and remarkably similar mass profiles. This is 
because magnetic field amplification happens after the 
disk is formed so the initial weak magnetic field will not 
affect significantly the dynamics of disk formation. In 
both cases, there are no sharp boundaries at the disk 
edge in surface density profile. From the mass profile 
one sees the disk radius to be ~ 3 kpc for both simula- 
tions. We overplotted the threshold gas surface density 
for st ar formation ^ 5 M^T^pc"^ derived from observa- 
tions (Martin & KennicuttJ |200il) • It is interesting to 
see that it is just near the disk radius. In the mass pro- 
file plot, the dark matter halo mass profile is also shown. 
It shows that the dark matter mass dominates baryon 
at r > 1 kpc. This explains why the rotational velocity 
profile is roughly fiat at r > 1 kpc with V(f) ^ 50 km/s, 
appropriate for the chosen NFW mass profile. In the ro- 
tation velocity plot, there is a big trough at ~ 1 kpc in 
both the Hydro and MHD curves. This is because of the 
fact that there are some big clumps at this radius whose 
spin provides a large cancellation to the mass-averaged 
rotation velocity. This also gives the trough in Toomre 
Q parameter discussed below. 

The bottom left panel of Fig. [3] shows the velocity dis- 
persions calculated by a = V< v'^ >— < > where 
< • > represents mass- weighted azimuthal average. The 
velocity dispersion has been argued to remain roughly 
constant in self-regulated regions of disks (|S ilk 1997), 
but this is still controversial (F erguson et al. I ll998). We 
can see that across the disk the velocity dispersions 
have about an order of magnitude fiuctuations centered 
around 5 — 10 km/s. So taking i t to b e constant is only a 
rough approximation. However. ISilk I ([l997 ) argued that 
the hot phase is the key in creating the self-regulated 
stage and determining the global velocity dispersions. 
Since hot phase is completely missing in our simulations, 
the situation might be different once supernova feedback 
is included. 

Fig, m shows t he Toomre Q parameter Q = Vtcj / (ttGT^) 
(|Toomre 1 1 19641 ), where Q is the angular velocity, a is 
the velocity dispersion, S is the surface density. In the 
original Toomre's formula, Q is replaced by the epicyclic 
frequency k, = {rdVt^ / dr + 4^7^)^/^. We used because 
our disk is highly turbulent, different averaging meth- 
ods for calculating Vt would give different n. In practice 
n ^ holds to within a factor of two. As can be seen 
from Fig. [3J throughout the disk Q fiuctuates around 
1, implying that the disk is Toomre unstable. Such a 
roughly constant Q suggests that the gravity-driven tur- 
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Fig. 5. — Mass-weighted joint temperature-B PDF at t = 1.088 
Gyr. The total mass in the disk is ~ 10^ M©, as can be seen from 
the mass profile plot of Fig. [3] The binsize is 0.02 in log space. 
The plus signs are the averaged magnetic field strength in every 
temperature bins. 
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Fig. 6. Mass-weighted joint density-B PDF at t = 1.088 Gyr. 
The solid straight lines are p oc B'^ and p oc B^-^. The color scale 
and binsize are the same as that of Fig. [5] The plus signs are the 
averaged magnetic field strength in every density bins. 



bulen ce in the disk is in a quasi-stationar y state at this 
time (|Gammie I I2QQTI : IWada et al.1 120021 ). This is also 
supported by the quasi-stationary density PDF (see sec- 
tion 13. 3p . Beyond the disk radius, Q rises sharply to 
> 1. 
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Fig. 7. — Top: mass-weighted joint density- PDF at t = 1.088 
Gyr. Bottom: mass-weighted joint density-S^ PDF at t = 1.088 
Gyr. The color scale and binsize are the same as that of Fig. [5] 
The plus sizes are the averaged By(Bz) in every density bins. The 
solid straight lines are p oc B"^ and p oc B^-^ in the top panel and 
p oc B"^ and p oc B^-^ in the bottom panel. 

3.2. The Growth and Structure of the Magnetic Field 

Comparing the density projections of the Hydro and 
MHD cases in Fig. (U we can see that the low filaments 
are more coherent in the MHD case. Correspondingly, 
the MHD disk has less low density regions. By com- 
paring the velocity field to magnetic field in Fig. [21 
we can see that the large-scale velocity field and mag- 
netic field are aligned in many places of the disk. This 
may be a result of self-organization in MHD turbulence 
by dynamic alignment of velocity and magnetic fields 
iBiskamp 1993). This implies that flows along field lines 
may be important for the formation of clouds ^. 

The edge-on slice of magnetic field lines shows that at 
large scale the poloidal field looks like a split monopole. 
This is the result of an initial uniform magnetic field 
along the rotation axis being dragged in by the col- 
lapse. The physics of magnetic field amplification in 

^ We thank Ralph Pudritz for alerting us to this point. 
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Fig. 8. — Top: time evolution of total disk gas kinetic energy 
(squares), thermal energy (diamonds) and magnetic energy (aster- 
isks). Bottom: time evolution of mass- weighted average plasma 
beta for total gas (squares), cold gas (diamonds) and warm gas 
(asterisks). 

MRI, i.e., the stretching of magnetic field lines increases 
their energy density at the expense of differential rota- 
tion, must also work here. From the face-on slices in 
Fig. [21 it can be seen that in the horizontal direction 
the magnetic field amplification is restricted mainly to 
the disk. But from the edge-on slice, it is clear that the 
field amplification extends significantly above the disk 
plane. As a result, there are highly magnetized bub- 
bles with /3 ~ 10~^ around the galactic disk. Those 
low beta bubbles form because density decreases sharply 
above the disk plane but the magnetic field strength de- 
creases much slower. Similar phenomena are well-known 
in the context of black hole accretion disks (e.g. Stone 
& Pringle 2001). The existence of such low beta bub- 
bles implies that the magnetic pressure may significantly 
influence the dynamics above the disk plane. Thus halo 
fields built i n this way may be crucial for halo-disk matter 
inter c hange (|Mc Kee fc Ostriker'1977': 'Norman fc Ike uchil 
1 19891 : iGomez d e Castro & Pudritz 1992), in me rgers 
with other galaxies (jToomre fc Toomrel 119721: [White! 
[1978[ ). ram-pressu re strippinpi; (|G unn fc Gott[ Il972[ ). 
galaxy h arassment ([Moore et al [[i996[ ) and galaxv stran- 
gulation (iLars on et al.[[l980[ ). 

Fig. [5] shows the mass-weighted joint temperature-B 
probability distribution function (PDF). From it we can 
see that the cold gas has a typical magnetic field strength 
~ 1 — 10 juG. This is remarkably consistent with observa- 
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Fig. 9. — Mass- weighted average of the toroidal magnetic field 
strength at t = 272 Myr (thin solid line), 544 Myr (dotted line), 
680 Myr (dashed line), 816 Myr (dash-dotted line) and 1.088 Gyr 
(solid line) 



tions in the Milky Way (|Crutched ll999'). Since f3 ex 
this imphes that in cold gas, (3 fluctuates by two order 
of magnitudes. Thus the average /3 gives only limited 
information for the magnetization of cold gas. Similar 
result has also b een found in local MH D simulations of 
molecular cloud ('Tillev fc Pu dritd [20071) . Furthermore, 
Fig. [5] also shows that magnetic fleld strength is similar 
in warm and cold medium (see also Fig. [2]). Then since 
thermal pressure is different in those two phases, this re- 
sults in the vast difference of plasma beta in those two 
phases. 

Fig. [6] shows the mass- weighted joint density-B PDF. 
It shows that the density-B relation roughly follows 
power law relation B ex rho^. For the high density disk 
gas, k = 0.5. This implies that the Alfven speed is sim- 
ilar in those density ranges. This is also consistent with 
the observ ations of B — p relation in dense molecular gas 
(|Crutcherlll999l ). The scatter in both observational data 
and our result are quite large. For the low density halo 
gas, k = 1.25. This results from the rise-up of the low 
beta bubble. 

To see how the magnetic fleld components scale with 
density. Fig. [71 shows the mass- weighted joint density-^^ 
and density- ^2 plots. It can be seen that both compo- 
nents have similar density-dependence. In the high den- 
sity gas, azimuthal magnetic component dominates over 
the vertical component, similar to what was f ound in 
shearing box simulations of MRI (jHawley et al. 111996) . 
In the low density bubbles, vertical magnetic fleld dom- 
inates over the azimuthal component. Interestingly, this 
is also consistent with radio polarization observations of 
nearby edge-on galaxies (Krause et al. 2006). 

To examine the ampliflcation history of the magnetic 
fleld, the top panel of Fig. [8] shows the time evolution 
of the total kinetic energy, thermal energy and magnetic 
energy for the disk gas. The magnetic energy increases 
exponentially by about flve orders of magnitude in the 
flrst 600 Myr. Afterwards it remains invariant with a sat- 
urated magnetic fleld strength in the disk. Some of the 
newly amplifled fleld after saturation will rise above the 
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Fig. 10. Density PDF for Hydro (top) and MHD (bottom) 
runs at t = 272 Myr (thin solid line), 544 Myr (dotted line), 680 
Myr (dashed line), 816 Myr (dash-dotted hne) and 1.088 Gyr (solid 
line). 

disk, magnetizing the halo and another fraction of it will 
be lost by dissipation. Furthermore, after saturation, the 
magnetic energy is still three orders of magnitude smaller 
than the total kinetic energy. This explains in part why 
in our simulation magnetic flelds did not influence signif- 
icantly the global disk properties. It could happen that 
at higher resolution magnetic flelds energy can be am- 
plifled to even larger value due to unresolved small-scale 
dynamo effects and may influence the disk more signifl- 
cantly. Nevertheless, it is encouraging that with the cur- 
rent resolution, the magnetic fleld strength in the cold gas 
already is very similar to the observed values in the Milky 
Way (Crutcher 1999). This is different from the result of 
dynamo simulations in MRI-driven turbulence where the 
total magnetic energy in the s aturated state is gen erally 
larger than the kinetic energy (jHawlev et al. mm . The 
primary differences are that turbulence in galactic disk 
is mainly driven by self-gravity. Furthermore, the simu- 
lated galactic disk has a multi-phase medium. The cold 
phase controls the saturation. And the warm phase still 
has a subdominant magnetic fleld in the saturated state. 
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This also contributes to the sub dominance of magnetic 
energy. 

The bottom panel of Fig. [8] shows the time evolution of 
the mass- weighted averaged logarithm of plasma beta for 
all gas, cold gas and warm gas in the disk. For t < 600 
Myr, the plasma beta in the cold gas quickly deceases 
from 10^ to one and then remains around unity. The 
magnetic field in the warm phase is always highly super- 
critical. This evolution behavior suggests that magnetic 
field amplification in disk galaxies is a self-regulated pro- 
cess. After the plasma beta reaches one in the cold gas, 
the mag netic fields be come dynamically important and 
buoyant (^ Parker 111 9661 ) . Further growth of magnetic field 
is quenched in the disk. This could be because of insuffi- 
cient resolution leading to too large dissipation, yet our 
52 pc resolution run gives very similar result. On the 
other hand, we can see clearly that magnetic fields begin 
to rise above the disk after saturation, forming the low 
beta bubbles. This suggests that at least up to 1 Gyr, 
buoyancy is a relevant mechanism of removing magnetic 
fiux from the disk. To further confirm this picture, we 
have run simulations that start with a ten times higher 
initial magnetic field, i.e. 10~^ G. They reach the self- 
regulating phase earlier and then drive more low beta 
bubbles magnetizing the halo in agreement with the pre- 
ceding interpretation. 

The thick line in Fig. [9] shows the mass-weigthed 
toroidal magnetic field at t = 1.088 Gyr. There are 
regular patterns of field reversal at late times. This is 
a consequence of field amplification by differentially ro- 
tating disks (Parker 1979). Fig. [9] also shows the time 
evolution of the mass-weigthed toroidal magnetic field. 
It can be seen that at early times the averaged toroidal 
fields are small, before they are amplified by disk rota- 
tion. The amplitude generally decreases with radius be- 
cause the rotation frequency decreases. A perhaps more 
surprising phenomenon occurred between t = 816 Myr 
to t = 1.088 Gry, when the toroidal field reversed almost 
exactly symmetrically. 

3.3. Statistical Properties of the ISM 

The density PDF is of interest with regards to star for- 
mation since star formation is thought to be de termined 
by the high-density tail of the PDF (lElmegree n 2002; 
iKrumholz & McKee "^OOS"; 'Wada & NormmJ 12007,) . A 
higher order (though not necessarily small!) effect is the 
spatial distribution of mass. When calculating the den- 
sity PDF, we considered only gas in the disk. Fig. [TOl 
shows the time evolution of this density PDF for both 
runs. It can be seen that the PDFs for the two cases are 
similar. The MHD calculation has slightly less volume 
in dense gas. In both cases the PDF reaches a quasi- 
stationary state within less than ~ 500 Myr. This quick 
establishment of a quasi-stationary PDF suggests that 
the gravity-driven turbulence in the disk reached a quasi- 
stationary and self-regulated state. This phenomenon, 
if robust even including supernova feedback, may play 
a crucial role in explaining the the observed universal 
relation between star formation rate and gas density 
(|Kennicutt I Il998l : IWong fc BlitTI l2002f ) . Our tempera- 
ture ffoor is a feedback on the ffow that in reality will be 
dominated by HII regions and supernova feedback. 

Our density PDFs cannot be fitted by log-normal 
distributions. At least, they do not show the typi- 
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Fig. 11. — Joint Temperature-Density PDF for Hydro (top) and 
MHD (bottom) runs at t = 1.088 Gyr. The color scale and binsize 
are the same as that of Fig. [5] The plus signs are the averaged 
temperature in every density bins. 

cal turn-over of a log-normal distribution. This is dif- 
ferent from previous simulations of disk galaxies which 
claimed that log-normal is a g ood fit to the densit y PDF 
(Tasker & Bryan 2006. [20071 : fWada Norman 1 [20071 ) . 
The crucial difference could be that those studies start 
from an isolated disk while in our case, the disk is built 
inside-out dynamically. Thus anytime in the evolution 
of the disk galaxy (< 1 Gyr) we have simulated, there is 
low density gas reaching the disk from outside. It is pos- 
sible that if we continue to evolve the disk long enough 
that after the initial halo gas all fell down to the disk it 
may also obtain a log- normal PDF. However, for a real- 
istic cosmological environment, infall should never stop 
completely. 

Fig. [11] shows the joint temperature-density PDF. 
Above p > 8 X 10~^^ g cm~^, the temperature is deter- 
mined by the fioor value which was introduced to ensure 
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Fig. 12. — Time evolution of the volume and mass fraction of cold (T < 10^ K) gas for the Hydro (thin line) and MHD (thick line) runs. 

the calculation resolves the Jeans length. So the tem- 
perature of the gas can at most be trusted below this 
density. From this figure, it is also clear that the ISM 
develops a two-phase structure with a warm T ~ 10^ K 
and a cold T ^ 300 K phase. However, unlike the origi- 
nally quasi-stationary picture of a two-phase ISM ( Field I 
[1965) , where the phase is characterized by a density and 
temperature determined by the balance of heating and 
cooling, in our case there are three order of magnitude 
density fluctuations in both cold and warm gas. This is 
due to the turbulent nature of heating in our simulation 
as discussed above. However, in both runs, one can still 
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cold phase. 

Fig. [12] shows the evolution of the volume and mass 
fractions of cold gas in both runs. It can be seen that the 
amount of cold gas is monotonically increasing in both 
volume and mass fractions. After 600 Myr, both runs 
reach a quasi-stationary volume fraction of ~ 0.8 and 
mass fraction ~ 0.99. 

3.4. Estimated Star Formation History 

Although we did not model star formation explic- 
itly, from the density PDF, we can estima t e the 
star formation rate (SFR) from (|Elmegreen~l l2002l : 
iKrumholz McKee II2005I : IWada Norman 11200 

M(> Per) 




600 
t (Myr) 



1000 



1200 



SFR = /, 



tdyn,t 



(2) 



where M(> per) is the total gas mass with density 
p > Per which can be computed from the density PDF, 
tdyn,cr = ^1 {GpcrY^^ IS the dynamical time at per and js 
is a constant specifying star formation efficiency (SFE). 

Using Per = 10-^^ g cm-^and fs = 0.01, Fig. [13] 
shows the star formation history (SFH) of the two simu- 
lations. It can be seen that while the SFR in the Hydro 
run is constantly increasing, in the MHD case, it reaches 
a peak at t ~ 800 Myr and then decreases. Since t ~ 800 
Myr is shortly after the time when the average beta in 



Fig. 13. — Star formation history of Hydro (thin line) and MHD 
run (thick line) with pcr = g cm~^and fs = 0.01. 

the cold gas reached unity, this suggests that magnetic 
forces provide further support in the cold gas after the 
ampliflcation saturated. 

Illustrating how sensitive the estimated SFR depends 
on per^ Fig. [T4l shows M(> pcr)/tdyn,cr as a function of 
Per- It can be seen that in both cases it does depend 
on Per- Thus to get the same SFR, the SFE fs should 
decrease in those densities ranges. However, it is inter- 
esting that the effect of magnetic ffeld is to make the 
dependence weaker. From n = 10 cm~^ to n = 10^ 
cm~^, in the MHD case M(> Pcr)/tdyn,cr increases by 
only a factor of 4. 

The current analysis of the SFR may be a poor pre- 
dictor of what happens when molecular cloud and star 
formation and feedback are included more realistically. 
However, the analysis in this section does suggest that 
after the magnetic growth saturates in the disk, magnetic 
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Fig. 14. — Ratio of M{p > pcr)/tdyn, cr a function of pcr 
computed from the density PDF for Hydro (thin hue) and MHD 
(sohd hue) runs at t = 1.088 Gyr. 

fields can significantly influence GMC formation which 
occurs in the strongly magnetized cold and dense gas. 

4. DISCUSSION 

Using the Spitzer values for viscosity and resistiv- 
ity, the magnetic Prandtl number of the warm ISM is 
~ 10^^. Thus, small-scale dynamos are expected to de- 
velop. However, those small-scale dynamo effects are not 
captured by our resolution. Never theless, previous sim- 
ulations (|Schekochihin et al. II2QQ4 ) have shown that for 
very high magnetic Prandtl number MHD turbulence, 
the growth of the small-scale magnetic fields is controlled 
by the rate of strain at the viscous scale and resulted in 
magnetic energy tending to pile up at the much smaller 
resistive scale. This is because for small scale dynamos, 
the magnetic field seems to retain its folded structure in 
saturation with direction reversals at the resistive scale 
(^Schekochihin et al. 2002). So small-scale dynamo pro- 
cesses may not be important for the large-scale growth 
of galactic magnetic field. At least, the situation is un- 
certain. Even if small-scale dynamos would contribute 
some fraction to the large-scale growth of the galactic 
magnetic field, our simulations of magnetic field growth 
would be a lower limit. Since we find that the field can 
be amplified to the observed value quickly, for the phe- 
nomenology of galaxy formation, it seems sufficient not 
resolve small-scale dynamos at the current level of de- 
tails^ 

Ide Avillez fc Breitschwerdt (2005) have discussed the 
results of 3D AMR MHD simulations of a local patch of 
ISM with sufficient vertical height to include the effect of 
halo-gas circulation. They find that gas transport to the 
halo is not significantly influenced by the magnetic field 
and they do not find the low beta bubble as seen in our 
simulations. We think the crucial difference of our sim- 
ulation from theirs is the presence of global differential 
rotation, which dominates magnetic field amplification 
in our simulations. 

An important limitation of the presented calculations 
is that we neglected supernova feedback while their dy- 
namics could be important for amplifying the magnetic 



field in the ISM (|Ferriere I fT992h . More importantly, 
super nova feedback will crea te the hot phase of the 
ISM (jMcKee Ostrike d flOTTh . which may play a cru- 
cial role in regulating star formation (|Silk I IiQQTI ). Fu- 
ture work needs to include this supernova feedback to 
study the importance of superbubbles on the structure 
of ISM, global star formation and magnetic field ampli- 
fication. Another limitation is that we neglected cos- 
mic ray pressure, which may be dynamically important. 
For example, cos mic ray pressure arising from the Alfven 
wave instability (|Kulsrud fc Pearce 1969) may prevent 
the downsliding of matter along the magnetic field line 
(|R,afikov fc Kulsrud l[2000h . 

5. CONCLUSIONS 

The simulations discussed in this work have shown that 
the amplification of magnetic fields and its large-scale 
topology is a natural part of the galaxy formation pro- 
cess. Our main findings can be summarized as: 

1. In the initial stage of the build-up of a galactic disk, 
the magnetic field is very weak and does not significantly 
affect the dynamics of disk formation. 

2. In a self-gravitating disk, vortex modes driven by 
the self-gravity of the gas grow rapidly and trigger the 
first generation of GMCs formation. 

3. In the cold and warm phase created by turbulent 
heating, the pressure is not isobaric. Temperature and 
density fluctuates by about three orders of magnitude in 
both phases. 

4. The dynamical formation of a galactic disk results 
in a global density PDF that cannot be fitted by a log- 
normal function. 

5. Differentially rotating galactic disks amplify a seed 
field of 10~^ G to microgauss levels in ~ 500 Myr. 

6. The growth of galactic magnetic fields is a self- 
regulated process. The saturation is reached first in the 
cold phase. After saturation, the field strength agrees 
with the observations of Milky Way magnetic fields. The 
cold phase, the magnetic field strength fluctuates by 
about an order of magnitude which gives larger than two 
order of magnitude fluctuations in the associated plasma 
beta. The saturated total magnetic field energy is three 
orders of magnitude smaller than the total kinetic energy 
in the case presented here. 

7. The saturation results in highly magnetized mate- 
rial around galactic disks at all but the earliest cosmic 
epochs. The halo magnetic fields may significantly in- 
fluence the halo-gas interactions, accretion of gas in disk 
galaxies, galaxy mergers and the interactions of galaxies 
with their environment. 

8. After saturation, the toroidal field in the disk dom- 
inates over vertical components while in the magnetized 
halo, vertical components dominate over toroidal com- 
ponents. 

9. Large scale magnetic field and velocity field are 
aligned at many places of the disk. This implies that 
cloud formation is likely channeled by flow along field 
lines. 

10. After saturation, the magnetic field strength is 
similar in the cold and warm medium. 

11. Global toroidal field reversals develops naturally 
in a differentially rotating disk. 

12. Magnetic fields can suppress star formation by pro- 
viding additional pressure support in the cold gas after 
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saturation. 

We have presented the first global numerical experi- 
ments studying magnetic field amplification and evolu- 
tion during the formation of disk galaxies. We find it 
encouraging that many of the aspects discussed here are 
consistent with evidence from observations (e.g. com- 
pare with Beck 2007). Being able to include more of 
the relevant physics such as molecular cooling, radia- 
tion and supernova feedback a realistic treatment of star 
formation, and a consistent treatment of cosmic ray ac- 
celeration transport, will render the next generation of 
numerical models useful for a more comprehensive un- 
derstanding of the physics of galaxies than the current 
state of the art. 
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